Real-life phylogeny analysis

In this notebook a single tree read from a newick string is analysed. We proceed by inferring the posterior via RjMCMC, followed by it’s analyses.

MCMC Controls

First we set up the MCMC parameters

set.seed(2)

burn_in <- 0.2
n_it <- 1e8
thinning <- n_it/1e4

Priors

Next, we parametrise and ready the priors for inference. We use the following parameters: - \(poisson(1)\) as the prior on the number of expansions - \(lognormal(3,3)\) as a relatively uninformative prior on parent population size - \(exp(10)\) as the prior on how long it takes for an expansion to reach half its carrying capacity in non-dimensional time units - \(lognormal(\log(N), 1/2)\) where \(N\) is the parent population size, as the prior on carrying capacities of individual expansions - \(gamma\) with mean \(\frac{N}{2}\) and variance \(\frac{N}{2}\) where \(N\) is the parent population size as the prior on when expansions happen N.B.: One should note that the model is conditional on the size of the parent population which effectively determines the tree scale. As such it is imperative that there is sufficient signal to infer this parameter with high enough confidence.

priors <- standard_priors(expansion_rate=1, 
                    N_mean_log=3, 
                    N_sd_log=3, 
                    t_mid_rate=5, 
                    K_sd_log=1, 
                    exp_time_nu=1/2, 
                    exp_time_kappa=1/2)

Run MCMC

We now read the tree and run the MCMC.

tree_in <- list.files(pattern=".*\\.nwk")
if (length(tree_in) > 0) {
  tree <- read.tree(file = tree_in)
} else {
  tree_in <- list.files(pattern=".*\\.tre")
  tree <- read.tree(file = tree_in)
}
tree <- makeNodeLabel(tree)
tree <- ladderize(tree, right = F)

if(length(tree$edge.length[tree$edge.length<0]) > 0) print("Negative branch lengths encountered and rounded to 1e-8")

tree$edge.length[tree$edge.length<0]<-1e-8

start <- proc.time()
expansions <- run_expansion_inference(tree, priors, 1, n_it=n_it, thinning=thinning)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%[1] "The acceptance rate is:  0.2108"
## 
  |                                                                            
  |=                                                                     |   2%[1] "The acceptance rate is:  0.2083"
## 
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%[1] "The acceptance rate is:  0.2067"
## 
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |===                                                                   |   4%[1] "The acceptance rate is:  0.2059"
## 
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%[1] "The acceptance rate is:  0.2052"
## 
  |                                                                            
  |====                                                                  |   6%[1] "The acceptance rate is:  0.2052"
## 
  |                                                                            
  |=====                                                                 |   6%
  |                                                                            
  |=====                                                                 |   7%[1] "The acceptance rate is:  0.2051"
## 
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%[1] "The acceptance rate is:  0.2048"
## 
  |                                                                            
  |======                                                                |   9%[1] "The acceptance rate is:  0.2046"
## 
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%[1] "The acceptance rate is:  0.205"
## 
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%[1] "The acceptance rate is:  0.2052"
## 
  |                                                                            
  |========                                                              |  12%[1] "The acceptance rate is:  0.2054"
## 
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%[1] "The acceptance rate is:  0.2058"
## 
  |                                                                            
  |=========                                                             |  14%
  |                                                                            
  |==========                                                            |  14%[1] "The acceptance rate is:  0.206"
## 
  |                                                                            
  |==========                                                            |  15%[1] "The acceptance rate is:  0.2064"
## 
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |============                                                          |  16%
  |                                                                            
  |============                                                          |  17%[1] "The acceptance rate is:  0.2064"
## 
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |=============                                                         |  19%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  21%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |===============                                                       |  22%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=================                                                     |  24%[1] "The acceptance rate is:  0.2064"
## 
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%[1] "The acceptance rate is:  0.2063"
## 
  |                                                                            
  |==================                                                    |  26%[1] "The acceptance rate is:  0.2063"
## 
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  27%[1] "The acceptance rate is:  0.2064"
## 
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |====================                                                  |  29%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |======================                                                |  31%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |======================                                                |  32%[1] "The acceptance rate is:  0.2067"
## 
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |=======================                                               |  34%
  |                                                                            
  |========================                                              |  34%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |========================                                              |  35%[1] "The acceptance rate is:  0.2064"
## 
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%[1] "The acceptance rate is:  0.2063"
## 
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%[1] "The acceptance rate is:  0.2063"
## 
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%[1] "The acceptance rate is:  0.2063"
## 
  |                                                                            
  |===========================                                           |  39%[1] "The acceptance rate is:  0.2064"
## 
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  40%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |=============================                                         |  42%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |==============================                                        |  44%
  |                                                                            
  |===============================                                       |  44%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |================================                                      |  46%[1] "The acceptance rate is:  0.2067"
## 
  |                                                                            
  |=================================                                     |  46%
  |                                                                            
  |=================================                                     |  47%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  48%[1] "The acceptance rate is:  0.2066"
## 
  |                                                                            
  |==================================                                    |  49%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  50%[1] "The acceptance rate is:  0.2065"
## 
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%[1] "The acceptance rate is:  0.2064"
## 
  |                                                                            
  |====================================                                  |  52%[1] "The acceptance rate is:  0.2063"
## 
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%[1] "The acceptance rate is:  0.2063"
## 
  |                                                                            
  |=====================================                                 |  54%
  |                                                                            
  |======================================                                |  54%[1] "The acceptance rate is:  0.2062"
## 
  |                                                                            
  |======================================                                |  55%[1] "The acceptance rate is:  0.2062"
## 
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%[1] "The acceptance rate is:  0.2061"
## 
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  57%[1] "The acceptance rate is:  0.2061"
## 
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%[1] "The acceptance rate is:  0.206"
## 
  |                                                                            
  |=========================================                             |  59%[1] "The acceptance rate is:  0.206"
## 
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%[1] "The acceptance rate is:  0.206"
## 
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%[1] "The acceptance rate is:  0.2059"
## 
  |                                                                            
  |===========================================                           |  62%[1] "The acceptance rate is:  0.2059"
## 
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |============================================                          |  63%[1] "The acceptance rate is:  0.2059"
## 
  |                                                                            
  |============================================                          |  64%
  |                                                                            
  |=============================================                         |  64%[1] "The acceptance rate is:  0.2058"
## 
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==============================================                        |  65%[1] "The acceptance rate is:  0.2058"
## 
  |                                                                            
  |==============================================                        |  66%[1] "The acceptance rate is:  0.2058"
## 
  |                                                                            
  |===============================================                       |  66%
  |                                                                            
  |===============================================                       |  67%[1] "The acceptance rate is:  0.2057"
## 
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  68%[1] "The acceptance rate is:  0.2057"
## 
  |                                                                            
  |================================================                      |  69%[1] "The acceptance rate is:  0.2057"
## 
  |                                                                            
  |=================================================                     |  69%
  |                                                                            
  |=================================================                     |  70%[1] "The acceptance rate is:  0.2058"
## 
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |==================================================                    |  71%[1] "The acceptance rate is:  0.2058"
## 
  |                                                                            
  |==================================================                    |  72%[1] "The acceptance rate is:  0.2058"
## 
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%[1] "The acceptance rate is:  0.2057"
## 
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  74%[1] "The acceptance rate is:  0.2057"
## 
  |                                                                            
  |====================================================                  |  75%[1] "The acceptance rate is:  0.2057"
## 
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  78%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |=======================================================               |  79%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  80%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |=========================================================             |  82%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  83%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  85%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |============================================================          |  86%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |=============================================================         |  86%
  |                                                                            
  |=============================================================         |  87%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |==============================================================        |  89%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  90%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  91%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |================================================================      |  92%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%[1] "The acceptance rate is:  0.2056"
## 
  |                                                                            
  |=================================================================     |  94%
  |                                                                            
  |==================================================================    |  94%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |==================================================================    |  95%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |====================================================================  |  97%[1] "The acceptance rate is:  0.2055"
## 
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%[1] "The acceptance rate is:  0.2054"
## 
  |                                                                            
  |===================================================================== |  99%[1] "The acceptance rate is:  0.2054"
## 
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%[1] "The acceptance rate is:  0.2054"
elapsed <-proc.time() - start

print(paste0(n_it," RjMCMC iterations completed. Time elapsed:"))
## [1] "1e+08 RjMCMC iterations completed. Time elapsed:"
print(elapsed)
##     user   system  elapsed 
## 186433.1      3.2 186463.8

Process MCMC data

    saveRDS(expansions, file = paste0(base_dir, "/expansions_s2.rds"))
    expansions <- discard_burn_in(expansions, proportion=burn_in)
  plot(expansions, mode="traces")

Inspect the base model parameters

We’re mostly interested in the model indicator trace and histogram, as this indicates the number of expansions. The population size trace acts as a good sanity check as well.

unique_dims <- unique(expansions$model_data$dim)
dim_counts <- sapply(unique_dims, function(i) length(which(expansions$model_data$dim==i)))
mode_dim <- unique_dims[which.max(dim_counts)]
if (mode_dim > 0) {
  plot(expansions, mode="summary", k_modes=mode_dim)
} else {
  plot(expansions, mode="summary")
}
## Warning: `mutate_()` was deprecated in dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
## Warning: Removed 130 rows containing missing values (geom_text).

## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange  gtable[layout]
## 2 2 (1-1,2-2) arrange  gtable[layout]
## 3 3 (2-2,1-1) arrange  gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[arrange]
print(paste0("The mode for the number of expansions is: ", mode_dim))
## [1] "The mode for the number of expansions is: 3"

Analyse the marginal corresponding to the number of expansion mode

We will look in detail at the marginal corresponding to the mode model, and inspect the top \(K\) most likely expansions

if(mode_dim == 0) {
  eval_events<-FALSE
  print(paste0("No events detected"))
} else {
  eval_events<-TRUE
}
plot(expansions, mode="modes", k_modes=mode_dim)